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1 Introduction 

In 1950 W. Kuhn and A. Katchalsky produced a fiber which contracts in response 
to changes in pH (Kuhn 1950). In their demonstration a single fiber immersed in a 
fiuid repeated raised and lowered a weight proportional to the hydrogen ion content 
of the solvent. It was proposed at the time that such a device could be used as a 
linear actuator. One of the major drawbacks however was the exceeding slow response 
time, on the order of minutes. Since then a number of innovations have made this 
artificial "muscle" worthy of further study. First, fibers have been produced which 
contract in seconds and even tenths of seconds (DeRossi 1987; Suzuki 1989). Second, 
some fibers support considerable loads, on the order of lOON/cm^ (DeRossi 1988). 
The contraction rates and forces have become comparable, and in some cases even 
exceeding, that of a human muscle. Third, the physics of fiber contraction (or gels 
as they are more generally called) have become well understood (Tanaka 1973a - 
87b). Finally, technological innovations such as robotics and implantable artificial 
biological organs have created a demand for such devices. A number of contractile 
gel devices have already been constructed, such as a robot gripper (Caldwell 1990), a 
multifingered hand (Toyota 1990), and an artificial urethral sphincter (Chiarelli 1988; 
DeRossi 1985; 86). 

The design of these devices from a practical engineering perspective, particularly 
with regard to dynamic modeling and control has only been considered recently (Gen- 
uini 1990). 

In this paper we will consider the design of a simple linear actuator based on 
contractile gel fibers as well as the design and simulation of single mechanical linkage 
controlled by two antagonist muscles. The first section will discuss the design of 
a simple linear actuator and the second will present a dynamic model along with 
approximate parameter estimates. The third section will introduce a nonlinear sliding 
mode controller which achieves desired trajectory tracking for model inputs. 




Figure 1-1. Two antagonist artificial muscles control a single link. 

2 Design 

Although there are many possible actuator designs, this paper will focus on the 
direct chemical to mechanical energy conversion through the control of the hydro- 
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gen ion content in the muscle fibers. The basic concept is to use high performance 
miniature valves to the control the inflow of acid or base to modulate pH and control 
contraction. A schematic diagram of the proposed actuator is shown in figure 2-f. 
Two input lines containing O.fM HCL and O.fM NaOH enter the muscle at the base. 
A miniature value placed into the base material controls the infiow of fiuid. An irri- 
gation system design to facilitate mixing is composed of numerous small tubes which 
are interspersed among the fibers. Finally, a drainage tube allows the waste fiuid (i.e. 
salt water) to be removed. 

The tubes and wetting surfaces of the valves are constructed from tefion which 
is chemically inert and the fibers are made from Poly- Vinyl Alcohol PVA. Although 
many different contractile materials exist (Tanaka 1990), PVA has been used suc- 
cessfully by a number of researchers and has a relatively high tensile strength. The 
tendons are made of Spectra and the tendon to PVA connections are machined Del- 
ran. The fibers are affixed at their terminations with epoxy which is also chemically 
inert. 
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Figure 2-1. Two fluid microvalues meter desired amounts of acid or base 
through an irrigation system to a bundle of contractile flbers. 

3 Model 

The dynamic model of the system is composed of three basic parts: the fiuid 
conveyance system, the hydrogel contractile fibers and the mechanical linkage. The 
following sections will address each of these systems in detail. 



3.1 Fluid System 

A schematic diagram of the fiuid system is shown in figure 3-1. It is assumed the 
inlet fiuid line is under a moderate pressure Po and a microvalve with controllable 
resistance Ry modulates the infiow of fiuid. The system of irrigation lines is modeled 
as a fiuid resistance Rt and an inertance /. Finally, the compliant sheath into which 
the fiuid fiows is modeled as a capacitance C and the exit line, a resistance Re. 
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Figure 3-1. The fluid model consists of a pressure source Pg at the inlet of 
the system. The valve, tubes and drain line are modeled as fluid resistances 
Ry , Rt and R^ respectively. Finally, a fluid inertance I is included to 
account for the mass of the fluid and a capacitance C to model the expansion 
of the actuator chamber. 



Valve 



Although there are commerciaUy available metering valves which continuously modu- 
late flow, they are generally too large from this application. Commercial piezo-electric 
or solenoid fluid values also exist in the appropriate size, although they have only two 
states, open and closed. It may be possible to use binary state valves in a pulse width 
modulation scheme to meter fluid, but this may introduce undesirable affects such as 
water hammer, excessive part wear and slower response. Alternatively there are some 
experimental metering valves and pumps which are of the correct dimension (in fact 
some are significantly smaller than this application demands creating the possibility 
of very small linear actuators). It may also be possible to use molecular valves in 
the form of biological or artificial membranes whose porosity is controlled by small 
voltages. 

In any case it is assume some mechanical valve of the appropriate dimension can 
be constructed, either a commercial two-state valve or an experimental multi-state 
valve. Based on commercial data sheets for fiuid resistance we can approximate the 
resistance of a valve by 

Ry = a{buy)''^, 

where Ry is fiuid resistance in N^/^sm""*, a = 4.09 X 10^ N^/^sm""*, 6 = 2.5 X 10"'* m/v, 
and n = —1.7925 are constants. For two-state valves we will assume Uy = Uo + sgn(M) 
and for multi-state valves Uy = Uo + u, where Uo is some small value and u is the 
control signal input, to 12 v. These parameters where roughly based on values 
given for the LEE Interface Eluidic 2-way (normally closed) microvalve. This valve 
is energized with 12 v at 250mW, switches in 1.5 msec and operates in a to 7 psi 
pressure range. The relation of pressure to fiow is given by 

where Qy is fiow rate in m"^/s, Py is pressure in N/m^ and Ry is fiuid resistance 

Nl/2gj^-4_ 



Irrigation System 

Multiple porous tubes attempt to uniformly enter and mix the fluid to increase the 

response time of the actuator. Consider a single tube as shown in figure 3-2. The 

tube has an inner radius r = 0.76 mm and a straight length / = fO cm after a 90° 

elbow of radius Te = 4r = 0.30 cm. Small holes of radius O.f mm are punctured in the 

tefion tubes. The holes are small enough so that an initial pressure P^ is necessary 

to induce fiow. This prevents back fiow into the acid-base lines, but for this initial 

model we will assume P^ is zero. The holes are place at 90° intervals about the cross 

section and every f mm along the length. 

Suppose the tube as we describe is non-porous and open at the end. If we assume 

a fully developed fiow under an entrance pressure of Po =1 kPa (above atmospheric), 

the fiuid fiow is 

ttP r'^ 

« = ^' '^) 

where fi = pp = 10""^ kg/m s (for water, 25°C), r = 7.60 X lO""* m, and / = O.I m. 
This yields a fiow of Q = 1.3 X I0~^ m"^/s and an average fiuid velocity of Uavg = 
Q/irr'^ = 0.72 m/s. The Reynolds number is Re = vd/p = 1089, where v = Uavg, 
d = 1.52 xIO""^ m and v = IO~^m~^/s, which is below the critical value Red,cr8t = 2300; 
therefore we may assume the fiow is laminar. For higher pressures, equation I yields 
Reynolds numbers in the turbulent range, but the significant resistance offered by the 
sealed porous tube reduce the fiow rates and maintains laminar conditions. As we 
will later show the maximum fiow rate is approximately Qmax = 1-65 X I0~^ m"^/s, 
which yields Uavg = 0.1 m/s and Re = 152, which is certainly laminar. Finally, the 
entrance region before the fully developed fiow is approximately /e = 0.06 Re d = lA 
cm; therefore we can assume the boundary layer is fully develop within the entire 
length of the tube. 

The wall shear force will result in a head loss and thereby create non-uniformity 
fiows through the small irrigation holes. However, since the fiow will decrease along 
the length due to seepage, the head loss will not be linear as in the case of homogeneous 
circular tubes. To approximate the head loss we will assume continuous porosity and 
a uniform fully developed fiow at each point x of tube. The wall shear stress is given 
by 

4:iJ.Q{x) 

which results in a pressure gradient 



Trr^ 



dP[x) —8fiQ{x] 
dx vrr'* 



If we assume a flow loss due to porosity 

dQ{x) 
dx 

we have the resulting differential equation 



dPix? 



-aPi 



„ - = hP{x), 
d/x 

where b = Sa/i/vrr'*. If we solve this differential equation with boundary conditions 
P(0) = Po and Q{1) = 0, we have the head loss as a function of length 
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2d 
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where c = v6 = 3.4 m ^. The head loss in the tube is not signiflcant in this case as 
shown in flgure 3-2. 
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Figure 3-2. Head loss due to wall shear force in the circular porous tube 
can generally be ignored. 

For the purposes of the dynamic model, it may be more more useful to repre- 
sent the entire porous tube system as a fluid resistance Rt- Following a simplifled 
convention, we will assume the fluid resistance is given by 

R = Ro + pi, 

where R is fluid resistance Ro is the resistance due to the entrance, p is the resistance 
per unit length and / is the length of the tube. The entrance resistance is given by 



R„ = ad 



(2) 



where d is the inner diameter of the tube in meters, a = 4.1 X 10"^ and n = —1.8. For 
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described above, the resistance per unit length is p = 5.3 X lO^N^'^sm"^, thus the 
total resistance for the tube is 5.13 X 10® N^'^sm"'*. Using equation ?? to compute 
the resistance of the irrigation holes, we have a resistance 4.0 X 10^° N^'^sm"'* for 
holes 0.1 mm in diameter. However with 80 holes/cm and 10 cm of tube, the total 
resistance is 1.6 X 10® N^/^ 



sm 
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For the acid and base systems we will have six irrigation tubes each, resulting in 
a total fluid resistance of Rj = 9.65 X 10^ N^'^sm"'* The relation between pressure 
and flow is again given by, 

RjQt = Pt- 

Inert ance 

Although the inertance of the fluid is small, fast switching times may cause its effects 
to be important. The inertance is given by 

,4' 

where p = 10"^ kg/m^ is the fluid density, / = 0.1 m is the length of the tubes and 
A = 1.1 X 10~^ m^ is the total cross sectional area. For this system, the fluid inertance 
is / = 9.2 X lO^kg/m'* and the dynamics are given by 

IQi = Pi. 

Suppose the fluid accelerates to 1.3 X 10~^ m"^/s in 10 ms under a pressure of 1 kPa. 
Then the inertance / would be 7.7 X 10^ kg/m'*. This is comparable to the value given 
for the model. So for switching times on the order of 1 ms, inertance effects become 
significant. 

Capacitance 

In this particular design, the fiber system is covered with a compliant sheath. Since it 
is fiexible, it can be approximated as a fiuid capacitance; that is, when fiuid is forced 
into the chamber the sheath will stretch and resist the fiow. If we model the sheath as 
a thin-walled uniform cylinder with closed ends, the elemental stress tensor is given 



by (Tr = 0, (Tg = Pr/t^ (Tz = Pr/2t and Tjj = 0, where P is the pressure, t is the wall 
thickness and r is the radius. Given the generalized stress-strain relations, 
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and the strain-displacement equations 



Ar = reg 
A/ = /e,, 



we can solve for the change in volume as a function of pressure 

V 2Et ; ' 

where E is the Elastic Modulus and p is Poisson's ratio. Therefore the fluid capaci- 
tance is 

7rr3/(3-4;/) 

~ 2Ei ■ 

If r = f cm, / = 10 cm, E = 10^ N/m^ (approximate for rubber), p = 0.5 and t = O.I 
mm, 

C = 1.6 X 10"^^ mVN. 

With a fluid flow of Q = 1.3 X I0~^ m"^/s, this capacitance will become significant 
if the pressure changes more than I MPa in a I ms. Therefore we will ignore the 
capacitance. 

Drain 

The drain line consists of a single tube de = 0.3 cm and a check valve. Employing 
equation ??, this results in a fiuid resistance of Rp = 2.22 X 10® N^'^sm"'*, for a drain 
line 20 cm in length. 

Fluid dynamic equations 



The elemental pressures are related by 

Po = Pv+Pi + Pt + Pe 

and the flows by 

Q = Qv = Qi = Qt = Qe- 

Substituting the elemental relations, we have 

Q = \[Po- [Rl + Rl + Rl)Q' 

where Q is the flow into the fibers. Note the steady state fiow rate for the open valve 
is 



P 



O 



Rl + Rl + Rt 

-7 ^3 , 



For a fully open value Qss = f -65 X fO m"^/s. Note the total volume of injected fiuid 



IS 



V = Q. 



Mixing and Diffusion 



When the fiuid (either acid or base) is injected into the fiber network there is a certain 
lag time before the entire system uniformly equilibrates to the new concentration. We 
could model this mixing as the diffusion of HCl or NaOH in water. If the tubes are 
uniformly spaced as in figure 3-3, we can draw concentric circles around each with 
radius Vo equal to half the distance between the tubes. The mixing time could be 
approximated, as the time for the fibers at radius Vo to experience a certain concen- 
tration of fiuid. As a simplification we can approximate the diffusion with the one 
dimension diffusion equations 



J 

dc 



or 



Q^2 



where J is the diffusion fiux per unit area, c is the concentration in moles per volume 
and D is a constant. The solution is the standard diffusion model 



1 — erf 



2VDt 



where 

2 
erf(x) = —= I e~^ dy. 



The boundary concentration Co^ of course, decreases as particles diffuse into the hquid. 
However, as a first order approximation to the diffusion time, we assume Co is constant. 
Then the concentration at Vo is one-half its initial value when t = lA ms for D = 
3.05 X fO-4 mVs (O.fM HCL at 25°C) and r^ = f.5 x fO'^ m. 




Figure 3-3. Arrangement of acid and base lines in the cross section of the 
muscle. 

The diffusion model, however, is highly inaccurate since the exit velocity from the 
pores in the tefion tubes is relatively high and together with the motion of the fibers 
will cause turbulent mixing, not diffusion. Therefore as a simple approximation, we 
will assume the fiuid mixing and diffusion is a first order lag, hence 

Tm(^ ~r C ^ Co, yo J 

where c is the concentration of acid or base in the solvent, r^ is a time constant 
which we will take to be r^ = f .4 ms (based roughly on the previous discussion) and 
Co = Vi/Vx ratio of solute to solvent, where Vt is the total volume of the actuator, 
Vt = 2.92 X fO"^ m^. 

3.2 Contractile Fibers 
pH 

The contractile fibers respond roughly linearly to pH, so it is important to calculate 
the pH in the solution based on the amount of acid or base entered into the chamber. 
An acid is a chemical species having the tendency to lose a proton and a base is a 
species tending to accept a proton. The quantitative measure of acid strength is the 
acid dissociation constant; that is the equilibrium constant for the reaction 



HA + H2O = H+0 + A" K 



+0 + A- A-S21I^ 
3" + ^ ^'- [HA] 



Along with this reaction there is always the self dissociation of water 

[H+0][0i7-] = /C 

If the original amount of acid A added to the solution is [HA]o, then the material 
balance relation is given by 

[HA]„=[HA]+[A-]. 

Also the solution must be electrically neutral, hence 

mO] = [A-] + [OH-]. 

Solving these equations for the hydrogen ion concentration yields the cubic 

[H3+0]([H3+0]-/C/[H3+0]) 



[HA]„-[H+0] + /C/[H3+0] 



K 



However, we may assume that the H||0 from water is negligible compared to that 
from the acid, hence [H||0] = [A~]. We may assume that for strong acids there is 
complete dissociation; that is [A~] = [HA]o, which yields the trivial result 

[HA]„=[H3+0]. 

Similarly, for base a B we have equilibrium 

B + H,0-HB+-OH A-='™™. 

mass balance, 

[B]„ = [B] + [+BH] 

and electric neutrality 

[BH+] + [H3+0] = [-0H]. 

Now we can assume total dissociation 

[B]o = [+BH], 

which is a valid except for very high concentrations. Also we can assume the contri- 
bution of ~0H from self-ionization of water is negligible, hence 

K K 

^ ^ ^ " [-0H] - [B]/ 
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Considering the strong acid and base together we can assume 

' [HA]„ - [B]„ [HA]„ - [B]„ > 



[H3+O] 



Kw / ( [B] o - [H A] ) otherwise 



Since pH is pH = — log[H;^0], we have 



pH 



-log([HA]„-[B]„) [HA]„-[B]„>0 
14 + log([B]o - [HA]) otherwise 



Diffusion 



Besides the lag due to mixing among the fibers, there is a lag due to the diffusion 
within the fibers themselves. This lag is due to bulk macromolecular rearrangement 
of the polymers in the network and is approximately proportional to the square of 
the cross sectional diameter. The time constant is 

T = cd^ , 

where d is the diameter in meters and c = 2 X 10^ s/m^. Thus for fibers 10 //m in 
diameter contract in approximately 0.1 s. This can also be approximated to the first 
order by 

rpH + pH = pHg, 

where r = 0.2, pH is the pH in the fiber and pH^ is that in the solvent. Notice this 
time constant is substantially slower than the one given in equation ??. Therefore 
we will ignore the lag due to mixing and consider only the diffusion within the fiber 
as the major delay. 

Hydrogel Physics 

A gel is a cross-linked network of polymers immersed in a fiuid, which can undergo 
a volume phase transition in response to changes in external conditions such as tem- 
perature, solvent, pH, light and electric fields. Three competing forces act on the 
polymer gel network: rubber elasticity, polymer-polymer affinity and hydrogen ion 
pressure. Competition between these forces, collectively called the osmotic pressure^ 
determine the equilibrium volume of the network. The phase transition has been an- 
alyzed in terms of the mean field theory of swelling equilibrium of gels. The osmotic 
pressure of a gel is given by Flory's equation. 



TT 



NkT 



+ ln(l 



+ 



2kf' 



+ vkT 




+ i^fkT ^ 
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where N is Avogadro's number, k is the Boltzmann constant, AF = AH — TAS 
represents the difference between the free energies of a polymer segment-segment and 
polymer-solvent interaction, and AH and AS are the enthalpy and entropy respec- 
tively. In equilibrium, the osmotic pressure is zero. Thus by solving the above equa- 
tion, we can theoretically compute the fiber force. However because of its complexity 
and because of the availability of empirical data, we can describe the contraction of 
the fiber with the simple model shown in figure 3-4. The intrinsic force per unit area 
is fm is approximately 

/„ = 1.0xlO^N/m'(7-pH^), 

assuming a total cross sectional area 2.92 X 10"'* m^, the intrinsic force is F^ = 29.2 
N(7-pHj). Also the intrinsic stiffness is 

k^ = 1.43 X lO^N/m. 

and viscous damping is roughly estimated as 



30Ns/m. 




^m 



Figure 3-4. Simple model of contractile fiber. 

3.3 Mechanical System 

A model of the mechanical system is shown in figure 3-5. Assuming linear com- 
ponents for this part of the model, we have 



J9 + {B + 2RHm)0 + 2R''kme = R{Fi - F2 



where 
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Figure 3-5. Single link controlled by two antagonist muscles. 

3.4 Summary 

Using the model as outlined above we have the following system equations 

K- = QtSgn{u) 
i>H, = (-pH,.-log(10-^ + K/VT)/T 

e = {-{B + 2RH^)e-2R^k^e + R{Fi-F2))/J. 



However before we design a controller for this system let us consider a number of 
simplifications. First, the lag due to the inertance of the fiuid is negligible compared 
with the dynamics of the rest of the system. This was verified in a full dynamic 
simulation of the system, although it is fairly obvious since the time constant of the 
fiuid fiow is on the order of 10 ms while the diffusion lag is about 100 ms. This 
assumption allows us to compute the fiow Q directly from the input u. Thus we 
can reduce the order of the system by two and take the input to be the fiow rate 
Qi. In order simplify this problem further let us assume a single-input single-output 
system, by taking the co-contraction of the actuators as an offset force and controlling 
position through feedback through a single actuator. With these initial assumptions 
let us write the system equations in the form 

X = f fx, u). 
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Specifically, 



where 



Xi = u 

X2 = aiX2 + a.ilog{a2 + 0.3X1) 

X3 ^ X4 

X4 = 04X4 + 05X3 + 06 + arX2 



Xi = y(cm ] 

X2 = pH 

X3 = 6 

X4 = 6 



anc 



ai 
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-33 s-i 


02 
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10-^ 


03 
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0.343 cm-3 


a4 
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-135 s-i 


as 
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-6457 s-2 


ae 


= 


2.89 X 10^ s-2 


ar 


= 


-3.87 X 103 s-2 



4 Control 

We wish to control the joint angle X3 = 6 using the simplified model presented 
above. We would also like the controller to maintain zero tracking error in the presence 
of modeling errors. Therefore, we will design a non-linear sliding mode controller. 

However before we attempt a construct a controller for this fourth order system, 
let us initial consider one further simplification. Based on the numerical estimates of 
the model parameters the damping ratio of the linear mechanical dynamics is ^ = 1.4. 
Therefore let us initially assume the inertia can be neglected, and produce a controller 
based on the reduce third order system, 

Xi = u 

X2 = aiX2 + a.ilog{a2 + 0.3X1) 

f 3 = 04X3 + O5 + OqX2 
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where 



X2 = pH 

X3 = 6 



and 



ai 
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02 
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Now the controllability canonical form, 

x('^) = /(x) + 6(x)u, 

where x = [a;,i, . . . ^x^'^'^'Y , can be derived by differentiating the equation for the 
joint angie three times. 

,2, 



X3 = a^aQ log(a2 -|- 033^1) -|- (oiOe + 0406)2^2 + 0.^X3 -\- 0504 



4 = aia6(ai -|- a4)iog(a2 -|- 033^1) -|- {aiaQ{ai -\- 04) -|- a^ae)x2 



n aiGe 

a^X3 -\- a^as -\- -u 

(a2 + 03X1) 

Let X = X — Xrf be the traclcing error on x and define the sliding surface to be 
S{t) = {x G 3J^|s(x,t) = 0}, where 

3(x,t)=(| + A) X, 

for some positive A. Hence 

s = X -\- 2\x -\- \ X. 

Thus maintaining zero tracking error is equivalent to the first-order stabilization prob- 
lem in s. In addition a bound on s implies a bound on the tracking error vector; that 
is if \s{t)\ < $, then |5;(t)(')| < (2A)'e. We must now design a controller such that for 
all X ^ S{t), 

2 j^'^^l^.O < v\-s{x,t)\, 
15 



for some positive rj. The estimation of / is given by / and the bounds on the error is 

1/ - /I < F. 
Similarly, the term b is estimated by 6 = V^min^max where < 6niin ^ b < 6max- If we 



let /3 = ^6max/&min, then 

and the control law 
for s given above, 
and 



u = b (u — ksgn{s)) 
k = /3{F + ,y) + {/3 - l)\ul 
u = —f -\- x^ — \x — \^x. 



In this case we chose A = 10 s~^ to be approximately equal to the natural frequency 
of the system while the estimate errors to be approximately 10% the nominal values. 

5 Simulation 

In order to get a sense of the model dynamics a IHz sinusoidal input was given to 
the system. A plot of the pH as a function of time is shown in figure 5-1 and the joint 
angle in figure 5-2. The sliding mode controller discussed in the previous section was 
implemented and simulated for a desired joint trajectory 

Xd = sin(27rt) -|- 1.5. 

Figure 5-3 shows the pH and figure 5-4 shows the desired and actual joint position. 
Notice desired tracking is achieved after approximately 0.5 seconds. Figure 5-5 shows 
the effect of neglected link inertia. Note the joint trajectory is offset by approximately 
10° . Fourth and fifth order sliding mode controllers were attempted to compensate 
for the link inertia and fiuid inertance, but were found to be numerically unstable 
in the calculation of the higher order derivatives. Numerical filters corrected this 
problem, though the simulation times became unacceptable long. A more judicious 
choice of numerical parameters could possibly correct this problem. 
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Figure 5-1. pH in the polymers fibers in response to a IHz input signal. 
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Figure 5-2. Joint angle. 
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Figure 5-3. Response of pH to control of joint trajectory. 
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Figure 5-4. Desired and actual (modeled) joint trajectories. Zero tracking 
error is achieved after one second. 
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Figure 5-5. Effect of neglected joint inertia. 

6 Conclusion 

A sliding mode controller produced reasonable tracking performance for input sine 
waves on the order of 1-2 Hz, but deviated significantly at higher frequencies. This 
seems reasonable based on the relatively long lag times introduced by the polymer 
diffusion. It has been possible more recently, to produce even smaller fibers (about 1 
//m), while using UV induced transectional cross-linking to increase strength. These 
and other innovations should allow faster response times for practical actuator devel- 
opment and controller design. 
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